% SimGreedModelLoopFinal.m  - many iterations of the numerical Greed model for generating Figure 2 and the Greed Model Part of Table 1
% Guardado and Pennings (2023) "Shock Persistence and the Study of Armed Conflict: Empirical Biases and Some Remedies" International Interactions
% Contact: Steven Pennings (spennings@worldbank.org; steven.pennings@gmail.com)
%       Replicates data for Figure 2; output Figure2data.xlsx, which used in ExcelGraphsTables.xlsx (which actually draws the graph)
%       Replicates Greed model data for Table 1; output Table1data_Greed.xlsx, which used in ExcelGraphsTables.xlsx (which actually produces the Table)
%       Calls solveGreedSS.m; GreedAR1.mod; GreedSV.mod
%       Uses Matlab R2018a (Optimization Toolbox; Statistics & Machine Learning toolbox); Dynare 5.3 (www.dynare.org)


%% Model set up %%
clear
clc
cd C:\Users\wb464833\Dropbox\SeasonalViolence\TheoryNote\Replication\GreedModel
addpath C:\dynare\5.3\matlab

% Parameters (we add the "1" as use this to keep variable later on")
global beta thetabar alpha phi psi  rho lambda 
beta=0.99;                  % Quarterly discount rate (in the manuscript we call this delta, because beta is the coefficient in a regression)
beta1=beta
thetabar=1;                 % Steady State TFP level(not logged!)
thetabar1=thetabar;
alpha=0.5;                  % Share of labor in production.
alpha1=alpha;
phi=3;                      % Opportunity cost is -phi (phi>1 to ensure a unique equilibrium)
phi1=phi;
psi=0.0065;                 % Scaling factor on the probability of winning the war
psi1=psi;
lambda=1;                   % Efficiency cost of government run by the rebels (1 is default; shouldn't be too low>0.5 or 0.6)
lambda1=lambda;
rho=0.966                   % Quarterly persistence of coffee prices (rho^4 annualized) - used in simulations in Figure 1A
%rho=0.983                  % Quarterly persistence of oil prices     (rho^4 annualized)
rho1=rho;

[vbar fval] = fsolve(@(x) solveGreedSS(x),[0.07]);
vbarstar=vbar;                                       % Steady state value of violence from Appendix 1.3.1
vbarstar1=vbar; 
p=psi*vbar^(1-1/phi)                                % Probability of winning 
pprime=(1-1/phi)*psi*vbar^(-1/phi)                  % Probability of winning derivative
wbar=alpha*thetabar*(1-vbar)^(alpha-1);             % 1. Definition of wages=MPL
vbarW=lambda*thetabar/(1-beta);                     % 3. Value Fn of winning  (lambda-=1)
vbarL=(wbar*(1-vbar)+beta*p*vbarW)/(1-beta*(1-p));  % 4. Value Fn of losing (out of power)  (lambda-=1)  
vbar1=vbar;
vbarW1=vbarW;
vbarL1=vbarL;

rng(0,'twister') % This is the default Dynare seed

%% Loop to generate Figure 2 in paper
%  First variation in the persistence of the shock 
rhoV_annual=[0;(0.3:0.01:0.86)';0.868851042654711;(0.87:0.01:0.92)';0.922262298269644; (0.925:0.001:0.99999)']
rhoV=rhoV_annual.^0.25;

N=length(rhoV); %141
reg_rhosimV=rhoV*0;
% reg_rhoIRFV=rhoV*0;   % regression on IRF generates same results
tic
for ii=1:1:N % takes 2 minutes to run
    ii
rho1=rhoV(ii);
dynare GreedAR1.mod noclearall
[B,BINT] = regress(v,[ones(1000,1) w]);
reg_rhosimV(ii)=B(2);
%[B,BINT] = regress(v_thetashock,[ones(40,1) w_thetashock]);
% reg_rhoIRFV(ii)=B(2);  % regression on IRF generates same results
end
toc

%  Then the seasonal shock
dynare GreedSV.mod noclearall
[B,BINT] = regress(v_thetashock,[ones(40,1) w_thetashock]);
B_SVregIRF1=B(2)
v_thetashock1=v_thetashock;
w_thetashock1=w_thetashock;

% Figure 2: Generating Figure2data.xlsx, which then ExcelGraphsTables.xlsx links to
Figure2data=[rhoV_annual, rhoV,reg_rhosimV, -phi1*ones(N,1), B_SVregIRF1(1)*ones(N,1)];
xlswrite('Figure2data.xlsx',Figure2data,'Figure2data')

%% Replicate  Table 1 (Greed model part; selected commodities) 
rhoV_annual_commodities=[0.863671000000000;0.956483000000000;0.899903000000000;0.868851000000000;0.949999000000000;0.970253000000000;0.922262000000000;0.958187000000000;0.967264000000000;0.969971000000000;];
N_commodities=length(rhoV_annual_commodities)
reg_rhosimV_commodities=nan(N_commodities,1);
for ii=1:1:N_commodities 
    ii
rho1=rhoV_annual_commodities(ii)^0.25;
dynare GreedAR1.mod noclearall
[B,BINT] = regress(v,[ones(1000,1) w]);
reg_rhosimV_commodities(ii)=B(2);
end

Table1data_Greed=[rhoV_annual_commodities, reg_rhosimV_commodities];
xlswrite('Table1data_Greed.xlsx',Table1data_Greed,'Table1data_Greed')
